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Abstract 

We introduce an inhomogeneously-nonlinear Schrodinger lattice, featuring a defocusing seg- 
ment, a focusing segment and a transitional interface between the two. We illustrate that such 
inhomogeneous settings present vastly different dynamical behavior than the one expected in 
their homogeneous counterparts in the vicinity of the interface. We analyze the relevant sta- 
tionary states, as well as their stability by means of perturbation theory and linear stability 
analysis. We find good agreement with the numerical findings in the vicinity of the anti- 
continuum limit. For larger values of the coupling, we follow the relevant branches numerically 
and show that they terminate at values of the coupling strength which are larger for more 
extended solutions. The dynamical development of relevant instabilities is also monitored in 
the case of unstable solutions. 



1 Introduction 

In the past two decades, the number of applications of discrete, nonlinear dynamical models has 
increased dramatically. A diverse set of applications has emerged, ranging from the nonlinear optics 
of guided waves in inhomogeneous optical structures 012] and photonic crystal lattices to 
atomic physics and the dynamics of Bose-Einstein condensate (BEC) droplets in periodic (optical 
lattice) potentials [5] |5] IE] and from condensed matter, in Josephson-junction ladders to 
biophysics, in various models of double-stranded DNA ^2 E21 This broad span of research areas 
and corresponding applications has now been summarized in a variety of reviews [131 1141 HTfl 1161 117) . 

A model that has drawn a particular focus among these areas of applications is the so-called 
discrete nonlinear Schrodinger equation (DNLS) |16j . This model was first proposed in the context 
of nonlinear optics, where it describes beam dynamics in coupled waveguide arrays |18l 119) . but is 
equally applicable in other settings, such as the dynamics of BECs confined in deep optical lattices 
|51 I2U|. The relevant model involves the nearest neighbor coupling between adjacent waveguides 
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(wells of the optical lattice in BECs) and the local nonlinear self-action induced by the Kerr effect 
in each waveguide (or the mean- field inter-atomic interaction in the condensate setting). 

Typically, the above setup is homogeneous in that all waveguides or wells are identical. However, 
recently there has been a surge of activity motivated by the experimental tunability of the proper- 
ties of individual waveguides/ wells. In particular, in the optical setting, the interaction of discrete 
solitary waves with structural defects was examined in |21| . while "non-standard" solitary waves 
(discrete gap solitons) were observed in binary waveguide arrays [221 [22] . This activity has been 
recently reviewed in |24j discussing various aspects of "optics in non-homogeneous waveguide ar- 
rays" . On the BEC side, there are also similar developments involving not only the (attractive or 
repulsive) localized "defect" action of a laser beam on the condensate 1231201) but also the potential 
of creating the so-called "superlattice" structures by means of the superposition of optical potentials 
of different periodicity [2T) . 

In this context of inhomogeneous nonlinear dynamical systems, we propose here a novel setting, 
which we illustrate to have drastically different dynamical behavior than that we would expect from 
its homogeneous counterparts. In particular, we impose a spatial pattern on the nonlinearity, hav- 
ing the form of an "interface" between a set of defocusing Kerr waveguides on the one end and a set 
of focusing Kerr waveguides on the other, separated by a "transient" layer (interface) of a waveg- 
uide bearing intermediate properties between the two segments above. This is, in some aspects, 
reminiscent to the recent proposition in the context of BECs of spatially dependent nonlinearities 
(see e.g. [211123 EH] and references therein). We show that this setting already presents a wealth 
of static and dynamical behavior which is very different than its homogeneous analog. 

We focus on the localized, solitary wave excitations in the vicinity of the interface. Starting from 
the so-called anti-continuum limit of zero coupling , we show that the existence and stability of 
the localized solutions in the vicinity of the interface can be quantified for low couplings by means 
of a perturbation theory using the coupling constant as small parameter. As the coupling between 
the sites near the interface increases, the phenomenology becomes drastically different, leading the 
relevant solution branches to a termination through saddle-node bifurcations that would be absent 
in the corresponding homogeneous limit. Perhaps equally surprisingly, the more extended multi- 
pulse solutions appear to survive for larger values of the coupling than the single (or smaller size) 
pulse waves, which is again contrary to what is expected from the homogeneous limit. In this 
stronger coupling regime, we investigate the properties of the relevant solutions through numerical 
bifurcation theory and linear stability analysis. We use direct numerical simulations to illustrate 
the manifestations of the dynamical instabilities of those among the solutions which are found to be 
dynamically unstable. We believe that this example illustrates the rich diversity of behavior that 
can be manifested in such inhomogeneous settings. 

Our presentation is structured as follows: in section 2, we present the setup and analytical results, 
while in section 3, we study the model numerically and compare with the analytical results. Finally, 
in section 4, we summarize our findings and present our conclusions. 
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2 Setup and Analytical Results 



We consider an inhomogeneous lattice model described by a discrete non-linear Schrodinger equation 
of the following form, 

i—7T = -CA 2 u n - d n \u n \ 2 u n , (2.1) 
at 

where C is the coupling between the adjacent sites of the lattice, A 2 u n = (u n +i + u n -± — 2u n ) 
is the discrete Laplacian. The evolution variable is z in the optical case and t in the BEC case; 
for notational simplicity, we use t in what follows. The nonlinearity coefficient d n (where n is the 
spatial index) is determined by the intensity-dependent refractive index (in the context of optics) 
or the s-wave scattering length (in the context of BECs) of each waveguide (or optical lattice well 
for BECs). The center site, do, is assumed to have an intermediate value slightly greater than zero. 
For all n < 0, the refractive index is set to the defocusing value d n = —0.9. For all n > 0, the 
refractive index is set to the focusing value d n — 1.1. Note that do is set to the average of these 
two values i.e., to 0.1; it should also be mentioned that the results reported below were found to be 
typical of similar choices of the d n profile. 

We focus our attention on standing wave solutions of the form of u n = e lAt v n , where A is the 
propagation constant in optics or the chemical potential in BECs, and v n is the spatial (time- 
independent) profile, satisfying the steady state equation: 

G(v n , C) = Av n - CA 2 v n - d n \v n \ 2 v n = 0. (2.2) 

It can be easily seen (see e.g. Refs. [321|33]) that, without loss of generality, we can restrict ourselves 
to the class of real solutions of Eq. I|2.2|) . In the anti-continuum (AC) limit, i.e., for C = 0, the 
solutions are immediately obtainable in the form: — {0, -j-}, provided that d n > for all n. 
Using this solution and setting to non-zero values only specific individual sites to the right of the 
zeroth site, we prescribe the configurations (in the AC limit) for the various branches that will 
be subsequently examined analytically as well as numerically. The selected configurations are as 
follows; lower first branch |1>: single excited site at n = 0; upper first branch |l,e>: excited sites 
at n = and n = 1 in phase; lower second branch |1, — e>: excited sites n = and n = 1 but 
out of phase; upper second branch |1, — e, — e>: excited sites n = with positive sign and n = 1, 2 
with negative signs. Following the same pattern, the remaining branches are: lower third branch 
|1, — e, e>; upper third branch |1, — e, e, e>; lower fourth branch |1, — e, e, — e>; upper fourth branch 
|1, — e, e, — e, — e>; lower fifth branch |1, — e, e, — e, e>; upper fifth branch |1, — e, e, — e, e, e>. 

On the numerical side, we analyse these branches using the pseudo-arclength continuation method 
34 . This allows us to trace the branches past fold points. In particular, given a solution (vq, Co) of 
the equation (|2.2|l G(v, c) = and a direction vector (vq, Co), one can determine (vi, Ci) by solving 
the following system of equations: 

G(v 1 ,C 1 ) = 0, 
(vt - voH + (Cr - Co)C — As = 0, 

where As is a (small) arclength parameter. We use Newton's method to solve the system in Eq. 
lt2~3l for (wl,Ci): 




G{v x ,Cx) 
{vi - voH + (Ci - C )Co - As 



(2.4) 
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The next (normalized) direction vector, (v\,Ci), can be computed by solving: 

(2.5) 




To examine the linear stability of the stationary solutions obtained as described above, we use the 
perturbation ansatz 

u n = e lAt (v n + ea n e-^ + 6„e^* 4 ), (2.6) 

where e is a formal small parameter. By substituting Eq. p.fijl into Eq. H2.1[l and dropping higher 
order terms, the following system of linear stability equations is obtained: 

uja n = -cA 2 a n + Aa n - 2d n \v n \ 2 a n - d n v 2 n b* n , _ ^ 

uj*b n = cA 2 b n - Ab n + 2d n \v n \ 2 b n + d n v 2 n a* n . 

The numerical solution of the ensuing matrix eigenvalue problem for the eigenfrequencies u> and 
eigenvectors {a n , 6* } can be then used to characterize the linear stability (more precisely the spectral 
stability) of the solutions. Since the eigenvalues (eigenfrequencies) of the underlying Hamiltonian 
system appear in quartets, to ensure a spectral instability it suffices for the above linear system 
to possess an eigenfrequency with a non-zero imaginary part. When the solutions are found to be 
unstable, we use a fourth-order, direct integration scheme to examine the dynamical evolution of 
the instability. 

Having presented the main framework and numerical methods, we now turn to some analytical 
results. Our analysis will be based on perturbation theory from the anti-continuum limit, using the 
coupling strength C as the small parameter. In particular, we expand the solution as: 

v n =v^ +CvM +0(C 2 ). (2.8) 

It is easy to check that the stability problem of Eq. (|2.7|) can be rewritten for the eigenvalues A = iuj 
in the Hamiltonian form 

JHip = Xtp, (2.9) 

where tp is the infinite-dimensional eigenvector, consisting of 2-blocks of {u n ,w n ) T (the superscript 
T denotes transpose), where a n = u n + iw n , b n = u n — iw n for the eigenvector equations 12. 7|) . J 
is the infinite-dimensional skew-symmetric matrix, which consists of 2-by-2 blocks of 

\Jn,m — 

and TL is the infinite-dimensional symmetric matrix, which consists of 2-by-2 blocks of 

o A 





7~tn.,,, , 

)n,' 

The matrices (C+) n , m and {C-) n ,m are, in turn, defined as: 

( £ +kn = 1 -3dnl£, (£-)„,„ = (£±) n>n+ l = (£±) n+hTl = -C. 

Similarly to the solution itself, the matrix 7i in the neighborhood of the AC limit, can be expanded 
as 

oo 

n = n {0) + ^2c k n (k} , (2.10) 

fe=i 
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where HS ^ is diagonal with two blocks: 

= ( ~ 2 ° ) , n E S, 7M = T " |> n£ AS, (2.H) 





where S denotes the set of excited sites. Notice that in the C — limit, each excited site corresponds 
to a pair of zero eigenvalues in equation (|2.9() , while each zero-site corresponds to a pair of eigenvalues 
at ±1. 

Choosing for simplicity of exposition (and without loss of generality) A + 2C = 1 in Eq. I|2.2[l . the 
solution of the leading perturbation problem in l|2.8() is governed by the following equation: 

[l-3d n ( W (°)) 2 ]«W=^ 1 +4°) 1 . (2.12) 

One can apply this, e.g., for the 2-site solutions such as |1, e > and |1, — e > , to obtain the leading 
order corrections: 



1 , fl 



2' V di 



»W (±J^), (2 .i3) 



,(!) 



= -2( ± V* ) ' (2 ' 14) 



where the sign inside the parenthesis corresponds to the sign of excitation of the site indexed inside 
the square root. Similarly, for 3 excited sites the expressions become 



„(!) 



4<w (2 - 15) 

-^ ± v* ± v* ) - (2 - i6) 

-i(±\/^)- ( 2 - 17 ) 



2 v v dr 

One can correspondingly generalize these expressions for an arbitrary number of excited sites. 

We now turn to the perturbed stability problem. The small perturbation of size C cannot ren- 
der the eigenvalues of order O(l) unstable. Instead, the potentially "dangerous" eigenvalues for 
instability purposes are those which are located at the origin of the spectral plane in the AC limit 
(corresponding to the excited sites, as discussed above). The perturbed form Hi of the matrix 
relevant to the stability problem can be easily seen (from the perturbative expansion) to assume 
the form 




W« = -ZdnWW ( _ ) , <U = H^ Un = -[;"], (2.18) 



1 
v 1 

while all other blocks of TLn}m are zero. If we consider the (linearly independent) eigenvectors 
corresponding to zero eigenvalues of H.q, f„, then it was proved in |33j that in order to obtain the 
leading correction to the (zero) eigenvalues of the original problem, it is sufficient to consider the 
reduced problem 

7WiC = 7 iC, (2.19) 

where 

(Mi) m , n = (f m ,H^f n ) (2.20) 
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is an N x N matrix, with N being the number of excited sites. Once the eigenvalues 71 of this 
reduced problem are obtained, then the perturbed eigenvalues of the original problem are given by 
the form A = VCAi + 0(C), where Ai = 

One can then directly compute the matrix A4, e.g. in the 2-site and 3-site cases, as well as more 
generally, to have the forms 

(*0„,« = -2d n v^v£\ {M) n;n _ x = -cob(0„_i - 6 n ), (M) ntn+1 = -cos(0 n+1 - 6 n ). 

One can then obtain the following predictions for the leading order eigenvalues of some of the 
branches discussed above (we only explicitly present these predictions in the 2- and 3-site cases) 

A = ±2.69003VC, (2.21) 
A = ±2.69003iVC, (2.22) 

for the (unstable) |l,e > and (stable, at least for small C) |1, — e > modes, respectively. In the 
3-site, we have: 

• For the branch |1, e, e >, two real eigenvalue pairs: 

Ai = ±1.2777UVC, A 2 = ±3.098986VC; (2.23) 

• For the branch |1, — e, e >, the same eigenvalue pairs as above but multiplied by i (hence the 
branch is marginally stable for small C); 

• For the branch |1, e, — e >, one real and one imaginary pair in the form: 

Ai = ±1.63075iVC, A 2 = ±2.428092^(7; (2.24) 

• Similarly, the branch |1, — e, — e > has the same eigenvalues as |1, e, — e > but multiplied by i 
(so it is also always linearly unstable). 

Notice that one can, in principle, expand this type of analysis to any other configuration of interest. 
We now turn to numerical results in order to examine the validity of these theoretical predictions. 

3 Numerical Results 

Figure ^ summarizes our essential numerical results, presenting the squared I 2 norm of the solution 
(physically, the power in optics or the rescaled number of atoms in BEC) P — J2 n \ u n\ 2 for the 
various branches that we examined in our computations. [For the explanation of the branches that 
are shown, the reader is referred to section 2]. There is a number of features in this bifurcation 
diagram which are in extreme contrast with the corresponding homogeneous limit of this system. 
Firstly, the single pulse branch in the vicinity of the interface already terminates for quite small 
values of C; in fact, it is the first branch to terminate in a saddle-node bifurcation with the two-site 
mode 1 1 , e >. This is the analog of what would be termed "the Page mode" in the setting of intrinsic 
localized modes (ILMs). As the coupling increases, the site with index n = 1 starts becoming excited 
for the single pulse branch eventually colliding (in configuration space) with the 2-site mode and 
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annihilating each other. In the homogeneous limit of a focusing medium both of these branches 
would survive for any C, up to the continuum limit of C — > oo. A similar phenomenology emerges 
for the so-called twisted mode branch of |1, — e > that, in turn, is also linearly stable for small C; 
for larger C it eventually collides with the branch |1,— e, — e > and disappears in a saddle-node 
bifurcation. The same is also true for the pair of |1,— e, e > and |1, — e, e, e > and for that of 
\l,—e,e,—e > and |1,— e,e,— e, — e > also shown in the figure. Another interesting general trend 
illustrated in this diagram is that the more extended the branch (i.e., the more sites participating 
in the nonlinear wave), the larger the coupling strength for which it persists. This is also contrary 
to what one would expect from the homogeneous limit, where multi-site solutions can only be 
continued to a finite coupling which is typically larger for more localized structures. 

Figure [21 illustrates the details of the lower pair of branches in Fig. ^ In particular, the top left 
panel shows the profile of the modes and their corresponding stability for the stable |1 > and 
unstable |1, e > solutions. The continuation of the branches (up to C = 0.15 where they collide and 
disappear) is shown in detail. The instability of the unstable two-site solution is investigated in the 
right panel through a direct simulation showing its breathing evolution. Panel (c) reports the result 
of the full numerical simulation (solid line) versus the theoretical prediction (dashed line) , both for 
the instability eigenvalue of |l,e >, and the profile correction imposed by Eqs. (|2.13(1 - (|2.14|1 . It is 
readily observed that, for small values of C, the agreement between the analytical predictions and 
the numerical results is very good. Of course, for larger C, the analytical results are expected to 
be less successful due to the significance of higher-order corrections neglected in our analysis. 

Figure|21is similar to Fig. |2 but for the second pair of solutions in Fig. ^ namely for |1, — e > and 
its corresponding unstable companion |1, — e, — e >. These branches disappear together in a saddle- 
node bifurcation for C — 0.575. While |1, — e, — e > is always unstable due to a real eigenvalue for 
any C, the twisted mode is stable for small C, but becomes unstable for C > 0.095 due to the 
collision of two imaginary eigenvalue pairs with opposite Krein signature, leading to a quartet of 
eigenvalues through a Hamiltonian Hopf bifurcation . The top left panel of the figure illustrates 
the solution profiles, the bifurcation diagram, and the results of the stability analysis. The top 
right highlights the breathing evolution of the twisted mode state. The two bottom panels compare 
the eigenvalue prediction of Eq. (|2.22|) with the full numerical result (dashed vs. solid lines) and 
the leading order correction for the two side mode, in the case of the left panel. In the right panel, 
the eigenvalue predictions (one real and one imaginary) for the \l,—e,—e > branch (dashed line) 
are also compared to the corresponding numerical results (solid line), obtaining once again good 
agreement. 

In Figure 0] we examine the third pair of branches of Fig. ^ namely (the stable for small C) 
|1,— e,e > and (the always unstable) |1,— e,e, e >. These branches, in turn, collide and disappear 
through a saddle- node for C — 0.725. The top right panel shows the prediction (dashed line) 
versus the numerical results (solid line) for the leading order eigenvalues of the |1, — e, e > solution 
(discussed in the previous section). The theory correctly captures, at small C, the existence of 
two imaginary eigenvalues and their C 1 / 2 bifurcation from 0, but is somewhat less satisfactory 
quantitatively in this case. This branch becomes unstable around C — 0.08 due to the collision of 
one of these eigenvalue pairs with one of opposite Krein signature, leading once again to a quartet 
of eigenvalues. The details of the subsequent dependence of this unstable eigenvalue on C, both 
for this branch and for 1 1, — e >, depend also on this size of the domain for reasons similar to 
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those discussed in jSSI - The two bottom panels show the unstable evolution of the corresponding- 
solutions, illustrating an interesting phenomenon particularly in the case of the |1, — e, e > branch. 
The dynamical evolution favors the tunneling of the excitation from the position of the interface to 
a nearby site (in this case, mainly to n — 2). That is, the interface displaces the solution towards 
a position where the environment is more conducive (being surrounded by focusing sites) to the 
existence of a localized pulse solution. 

4 Conclusions 

In the present paper, we have introduced a new setting for the study of the recent theme of inho- 
mogeneous nonlinear lattices in nonlinear optics (waveguide arrays) and Bose-Einstein condensates 
(optical lattices and superlattices) . The setting consists of an interface between defocusing and 
focusing (repulsive interaction and attractive interaction, respectively) regions and a transient layer 
between the two. 

We have focused specifically on the statics and dynamics of coherent waveforms in the vicinity of 
this interface, and we have found that their properties are dramatically modified in comparison 
with those expected from the homogeneous case. Some manifestations of these differences can 
be quantified in the termination of the principal pulse branch (for small couplings) or the more 
prolonged (in parameter space) persistence of more extended structures in comparison with more 
localized ones. Furthermore, we have shown that the interface may induce a dynamical tunneling 
of the structures towards locations more favorable for their existence. We have also developed a 
systematic methodology based on the adaptation of the considerations of j^S] to inhomogeneous 
settings and illustrated how to use these to develop a perturbative treatment of the problem with 
excellent qualitative and good quantitative agreement with the full numerical results. 

There are many interesting questions that are suggested for the interface problem we have in- 
troduced. A prominent one concerns the dynamical evolution of localized structures towards the 
interface and their interaction (transmission, reflection or trapping) with that region. 
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0.1 0.2 0.3 0.4 0.5 0.6 0.7 

Figure 1: Bifurcation diagram of the first five branches. Plot of the solution's power vs. the 
continuation parameter, C. Dotted lines represent unstable regions. Solid lines represent initially 
stable regions. 
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Figure 2: (a) Profiles of wave configurations and eigenfrequencies at the value of C = 0.1 where 
the solid vertical line crosses branch one in the center diagram. Upper left: wave configuration 
from lower half of branch one, |1 >. Lower left: eigenfrequencies of the linearization around this 
solution. Upper right: Wave configuration from upper half of branch one, |1, e >. Lower right: 
eigenfrequencies from the corresponding linearization, (b) Top: evolution of max square modulus 
of the solution taken from upper branch at C = 0.1. Middle: spatial profile of the square modulus 
of u taken for c = 0.1, after propagation by 250 units. Bottom: space time contour plot of the 
square modulus of the solution, (c) The top panel shows the most unstable eigenvalue of the two-site 
solution of |1, e > as a function of the coupling strength C . The solid line is the full numerical result, 
while the dashed line denotes the analytical prediction. The bottom panel shows the correction for 
the central site n = and its neighboring site n = 1 given by the first order theory (dashed line) 
versus the corresponding numerical result (solid line). 
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Figure 3: (a) Similar to Figure 2 for branch two with profiles at C = 0.135, for the branches |1, — e > 
and |1, — e, — e >. (a)-(c) are similar as above, while panel (d) shows the eigenfrequencies (one real 
and one imaginary, as predicted by theory) of the mode |1, — e, — e > from the numerical results 
(solid) against the analytical predictions of section 2 (dashed lines). 
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Figure 4: Panel (a) shows the profiles and eigenfrequencies of the third branch (|1, — e, e > on the left 
and |1, — e, e, e > on the right) for C — 0.465. The right panel shows, for the 3 site mode |1, — e, e >, 
the dependence of its two eigenfrequencies against the analytical predictions (dashed lines). The 
bottom panels show, for each of the branches and for C = 0.465, the dynamical evolution of the 
principal sites, as well as the space-time contour plot of the square modulus. 
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